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We use computer simulations to study the onset of collective motion in systems of interacting 
active particles. Our model is a swarm of active Brownian particles with internal energy depot 
and interactions inspired by the dissipative particle dynamics method, imposing pairwise friction 
force on the nearest neighbors. We study orientational ordering in a 2D system as a function of 
energy influx rate and particle density. The model demonstrates a transition into the ordered state 
on increasing the particle density and increasing the input power. Although both the alignment 
mechanism and the character of individual motion in our model differ from those in the well-studied 
Vicsek model, it demonstrates identical statistical properties and phase behavior. 



I. INTRODUCTION 

Dynamic self-organisation and, in particular, mecha- 
nisms of swarming behavior of microorganisms, cells, and 
animals remain one of the most intriguing problems at 
the interface of physics and biology. Numerous phys- 
ical models of interacting self-propcllcd particles have 
been proposed recently to study these phenomena (sec 
review papers ) . All these models capture the essen- 
tial prerequisites for swarming: out-of-cquilibrium state, 
which is manifested in the self-propulsion of particles or 
other mechanisms of transforming external energy into 
directed motion, and aligning or attractive interactions 
between the particles. The motion of individuals has 
been described in the simplest case by particles moving 
with a constant speed and subjected to angular noise 0- 
. More advanced presentations of active agents can in- 
clude friction, thrust force, and noise, like, for example, 
the active Brownian particle model (ABP) or even 
a very detailed mechanics of cell or animal locomotion 
@j flo| . The interactions required for the transition from 
individual to collective dynamics have been introduced 
in a variety of ways. In the standard Vicsek model, the 
swarming results from the action of a collision-type in- 
teraction that aligns the velocity of each actively moving 
particle in a big ensemble to the average local velocity 
1^,0]. Alternatively, the particles' individual direction of 
motion can be coupled to the mean orientation or posi- 
tion of the swarm jll| - [l4| . In some implementations, the 
type of many-body interaction depends on the distance 
between neighbouring particles [l5|. 

The ABP model, based on the Langevin equation for 
particle velocity, has received much attention in literature 
and has been successfully applied to a variety of problems 
0. Depending on the type of system under study, dif- 
ferent types of coupling between the ABPs been used. 
Several realisations of the model assumed only conserva- 
tive interactions d, [l^ or chemical interactions [O, [3 
between moving agents. Another development was based 
on the theory of canonical-dissipative system [l^ . Erd- 
mann and Ebeling studied the active Brownian particle 



model with Oseen-typc hydrodynamic interactions [20[ 
and observed several swarming modes. Recently, Gross- 
mann et al. [2l| studied the onset of the collective motion 
in a system of ABPs with velocity alignment and both 
passive and active fluctuations and found not only ori- 
entational order-disorder transition but also bistable dy- 
namics states. Swimming particles with hydrodynamic 
interactions were studied using the mesoscale simulation 
methods or direct solution of the Stokes equation also 
showed an onset of collective dynamics [22, i23| . 

The special attention paid to hydrodynamic effects in 
these studies is not accidental: the language of hydro- 
dynamics is conceptually well suited for description of 
the swarm motion. The analogy has been explored al- 
ready in the early papers by Toner and Tu P, y, 113, HI] , 
who introduced a Navier-Stokes-like continuum model 
for active materials. The hydrodynamic behavior of ac- 
tive swarms can be inferred directly from microscopic 
description. Bertin et al. derived hydrodynamic equa- 
tions governing the density and velocity fields from the 
microscopic dynamics for a gas of self-propelled parti- 
cles with binary interactions [l^. One can notice that 
similar ideas are exploited in the mesoscale methods in 
fluid modelling such as multi-particle collision dynam- 
ics (MPCD) (23, Hill where a collision operator is used 
to align particles to the average local flow direction, or 
the dissipative particle dynamics (DPD) [1^, HO], where 
the hydrodynamics comes in through inelastic collisions 
between the particles. Both of these methods designed 
to respect the momentum transport (long- wave hydrody- 
namic modes) and to punish the fluctuations (the high 
frequency modes) to achieve the hydrodynamic behavior 
at longer time and lengthscales. Obviously, the swarming 
behavior can be achieved through MPCD or DPD-like in- 
teractions as well. In fact, the possibility to develop the 
collective dynamics through dissipative interactions has 
been recently investigated by Grossman et al. [3l|. In 
their model, the system of active particles with spring- 
dashpot interactions demonstrated a discontinuous tran- 
sition into the aligned state upon reduction of noise and 
various types of collective migration or vortex-like motion 
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depending on the confinement. Thus, we are convinced 
that a system with dissipative interactions can have at 
least quahtatively similar dynamics to the models with 
aligning interactions like the Vicsek model. In the spirit 
of these observations, it is tempting to test whether the 
quantitative features of the swarming behavior can be 
reproduced in a generic dissipative model. 

In this paper, we study the dynamics of such a model 
and demonstrate that collective motion regimes can be 
achieved in the same way as in the standard Vicsek 
model. We combine two well-developed approaches: the 
active Brownian particle model, which allows us to in- 
troduce the self-propulsion and interactions with the en- 
vironment in a transparent way, and dissipative interac- 
tions for the active particles, so that the collective dy- 
namics would arise from explicit pairwise forces. By an- 
alyzing the statistical properties of this hybrid model, we 
show that it has the same universal properties across the 
order-disorder transition as those reported for the Vicsek 
model despite the big differences in the implementation. 
In Section HI] we describe the construction of the model, 
in Section [ml we show its statistical properties in a wide 
range of parameters, calculate the phase boundary for 
the orientational order-disorder transition. We discuss 
the results in Section |IVl and then conclude the paper in 
Section El 



II. MODEL AND SIMULATION SETTINGS 

To study the dynamic self-organisation of active par- 
ticles we introduce a two level model. At the single par- 
ticle level we include the factors determining the particle 
motion in a thermal bath (such as a water reservoir): 
bath temperature and thermal noise/fluctuations, fric- 
tion, and a motor. At the higher level, we introduce dis- 
sipative interparticle interactions. At the two-body level, 
we also include some noise, whose nature is, however, dif- 
ferent from that of the environment. The noise at the 
many-body level refers to biomimetic behavioral features 
like imperfect alignment of particles to their neighbours 
or sudden voluntary switching of the direction of motion. 
Although the noise in the interaction is of non-thermal 
origin, it can also be characterized by some effective tem- 
perature. We can show that the collective behavior can 
be associated with these temperatures. 



For simplicity, we will always set the particle mass to 
unity. The first two terms in Eq. ([l} are standard 
Langcvin friction and random force, where 7^ is the co- 
efficient of viscous friction, which is set by the properties 
of the environment. The last term is the random force of 
strength and ^{t) is representing a Gaussian white 
noise with zero-mean and unit variance. The strength of 
the noise is set by the fiuctuation-dissipation relation at 
the ambient temperature 



(2) 



We should stress that the fiuctuations introduced in Eq. 
P]) do not depend on the particle's speed or direction of 
motion and therefore can be referred to as passive [2l| . 
The thrust term in the depot model has the form: 



dq 



dV/ 



(3) 



where d is the constant determining the rate of conver- 
sion of internal energy of the active agent into kinetic 
energy, c is the parameter controlling energy loss, and q 
is the constant determining the gain of energy from the 
environment. The steady state motion of the active par- 
ticles is characterized by the velocity Vq ~ Vq^ + V^y, 
which is expressed through the system's parameters as 
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c 
d 
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a,t q > ^^c/d [1, The steady state velocity distribu- 
tion for various q is shown in Fig. [TJ At (j = we observe 
the Maxwell distribution of the velocities corresponding 
to the systems's temperature , while at non-zero en- 
ergy influx rates we sec either a broadened distribution 
centered at zero (at q < -f^c/d) or two bell-like peaks 
around the stationary velocity Vqx = Vb/v^ with the 
peak width controlled by the temperature. Thus, the 
model exhibits a transition from dissipative to driven 
regime upon increase of the energy conversion coefficient 
q d, [1^, [i^l . In the driven regime, attractive and aligning 
interparticle interactions can lead to collective dynamics. 



B. Interparticle interactions 



A. Equation of motion of a single self-propelled 
particle 

At the single particle level, we follow the ABP model 
Q. Here, we consider the motion in two dimensions. 
The equation of motion of an individual particle i is the 
Langcvin equation for the velocity with an added thrust 
term 

= -7^V, + Ff + (1) 



The collective behavior is impossible in a system of 
non-interacting agents. While one can expect some 
swarming (particle clustering) already with isotropic cen- 
tral interactions, the global symmetry breaking and on- 
set of directed transport requires that particle velocities 
are aligned. Here we should note that a spontaneous 
transition into a globally aligned state is impossible in 
equilibrium system with perfectly elastic collisions and 
without any dissipation due to conservation of the total 
linear momentum. The local alignment can be realised 
by different means, the best known example being the 
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FIG. 1: Instantaneous ID velocity distributions for the active 
Brownian particles at c = 0.8, d = 2.0, = 0.3. Inset: Bi- 
furcation diagram for the stationary velocity Vb as a function 
of q. 



Vicsek model, where particles are aligned to the local 
mean velocity field [sj. Another example of the align- 
ing interaction is the hydrodynamic interaction of fluid 
molecules, solute particles, or swimmers [2^. At the mi- 
croscopic level, the onset of hydrodynamic behavior is 
achieved by suppressing the relative motion of the neigh- 
bouring particles with a friction force and preserving the 
local mean velocity. As a result, the fluid quickly relaxes 
to the stationary state. This idea is realised in a num- 
ber of mesoscale simulation methods, which are known to 
produce correct hydrodynamics: lattice Boltzmann (LB) 
method m, MPCD, and DPD. In the LB and MPCD 
implementations the collisions are collective, similar to 
the Vicsek model, while in DPD the friction force is pair- 
wise and is applied to each pair of colliding particles [33| . 
So, the latter method is ideally suited for our purpose as 
it presents a simple way to control the strength of the 
aligning interaction. 

Here, we introduce a dissipative force between the 
ABPs in the same way as it is done in the DPD method. 
The total force 'Fi{t) acting on each particle is then given 
by: 



(5) 



where Ff is the force that comes from interactions within 
the swarm. Ff consists of three parts: 
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where F^^ ^ , and F^ represent the conservative, dis- 
sipative, and random forces between particles i and j, 
respectively. The conservative force that reflects the ex- 
cluded volume interactions is defines as: 



Ti — Tj is the distance between particles i and j, 
|r.y| is its magnitude, and = Yij/vij is the unit 



vector from j to i . 
repulsion: 



F^{r) 



We choose F^ 



a 1 - 



to describe a soft 



0, 



r <rr 
r > ?> 
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where a is a parameter determining the maximum repul- 
sion between the particles, = fc/2 is the radius of the 
repulsion zone, and Vc is the cut-off distance. 

The dissipative force punishes the velocity differences 
between the neighbouring particles and, therefore, pro- 
vides a mechanism of relaxation of the velocity field to- 
ward the stationary state. We take it in the form of a fric- 
tion force applied to the component of the motion in the 
direction of the particle connecting vector, i.e. a speed 
adjustment for particle moving together in the same di- 
rection. 
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Similarly, the friction can be applied to the motion per- 
pendicular to the connecting vector [s^ ] , in which case it 
will predominantly act as an aligning interaction. In both 
cases, the parameter 7"^ controls the dissipative strength 
of the interaction and by varying it we can accelerate or 
delay the alignment. 

In a thermalized system, the stochastic force F,^: com- 
pensates the loss of kinetic energy due to the dissipative 
force. It provides random "kicks" in the radial direction 



causing misalignment of particles' velocities. 
F«(t) = y2^«(r.,)e,,(t)f., 



(10) 



where determines the strength of stochastic interac- 
tions, and ^ij(t) is a random variable with a Gaussian 
distribution and unit variance. In hydrodynamic simula- 
tions, it is usually required that the noise ^ij is symmetric 
in ij , the kicks satisfy Newton's third law and conserve 
total momentum [35| . This requirement, however, can be 
omitted for active particles. A swarm temperature, , 
can be defined via the fiuctuation-dissipation relation 
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In Eqs. ©-([I 



(r) and uj^{r) are weight functions, 
which for simplicity can be selected to be similar in form 
to the repulsive scalar function F'~^{r): 



r < Tr 



(12) 



(7) 



where F'~^ (r) is a non- negative (repulsive) scalar function 
determining the distance dependence of the repulsion. 



In this model, we can regulate the interparticle interac- 
tion by changing the effective temperature, T'^, and the 
friction coefficient, 7'^. The effective temperature deter- 
mines the average degree of alignment the system can 
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tolerate, while the friction coefficient determines the dis- 
sipative power of a single collision and the speed of relax- 
ation toward the stationary state. Note that in this case 
the friction and the noise depend on the particle relative 
position and velocities. The noise is generated by the 
individual itself so the fluctuations can be referred to as 
active [2l|. Clearly, the global ordering should depend 
on both types of fluctuations, thermal (environmental) 
and non-thermal (collective). The whole set of the DPD 
terms reflects the behavioural contributions to the mo- 
tion. For an animal or robotic systems it amounts to re- 
specting the excluded volume and adjusting the motion 
to the neighbours. The stochastic term in this context 
plays a role of angular noise or imperfect alignment of 
the agents to their neighbours' direction of motion. 



C. Simulation settings and motion statistics 

We used a two-dimensional system with periodic 
boundary conditions. The primary box size was fixed 
at 200 X 200 units and we varied the number of particles 
in the main cell to set the number density p. Simulations 
were performed with time step of At = 0.01. Particles 
were propagated using the Verlet algorithm (36j : 

r,{t + At) = 2r,{t) - r,{t - At) + At^F,{t) (13) 

The velocities of particles were calculated using Stormer- 
Verlet method: 



v.(t) 



r,{t + At) -ri{t- At) 
2At 



(14) 



Total number of time steps in each run was set to 1 x 10^. 
All simulations were performed with the following set of 
key parameters: Vc = 3, 7^ = 1, T^' = 0, 7^ = 0.3, 
= 0.3, a = 1, d = 2, c = 0.8. 
To characterise the collective motion in our model 
we use three different correlation/distribution functions. 
The velocity autocorrelation function is calculated as 



V,(0)-V,(<) 



A^\^|V.(0)|.|V,(0|/ 
The radial distribution function is defined by 

N N 



(15) 



9(r) 



N{N 



-T)(EE'5(--|r«.i)) (16) 



where S is the Dirac delta function, (•) stands for the en- 
semble average. The two-point velocity correlation func- 
tion is calculated as 



N N 



C\\ir) 



v.{t)-v,{t) 

N{N^1) \V^{t)\■\VAt)\ 



V- 1) 



(17) 



where i and j label particles separated by distance r = 
\rij\. With this definition, two particles with parallel 



(antiparallel) velocities give a correlation of -1-1 (— 1). 
The angular brackets denote the ensemble average. To 
characterise the swarming behavior of the particles we 
also perform a cluster analysis. Cluster in our model is 
defined as a group of particles with a distance between 
neighbours smaller or equal to the cut-off radius there- 
fore, particles interacting directly or via neighbouring 
agents are included into one cluster. We calculate the 
number of clusters and mean cluster size. 

We characterise the orientational ordering by the po- 
lar order parameter, which quantifies the alignment of 
the particle motion to the average instantaneous velocity 
vector 



ifiit) = (cos 0,(0) 



1 ^ 
-T 

N ^ 



^^|V,(t)||(V(t))| 



(18) 



This order parameter has been extensively used to de- 
scribe the orientational ordering in various systems of 
self-propelled particles 0, [s^]. It turns zero in the 
isotropic phase and finite positive values in the ordered 
phase, which makes it easy to detect the transition. 

To locate transition points precisely we also calculated 
the Binder cumulant 1381 



Gr = 1 - 



(19) 



where stands for the time average and L denotes the 
value calculated in a system of size L. The most im- 
portant property of the Binder cumulant is a very weak 
dependence on the system size so takes a universal 
value at the critical point, which can be found as the in- 
tersection of all the curves Gl obtained at different sys- 
tem sizes L |37| at fixed density. To detect the transition 
points in g — p plane precisely we plot three curves for 
different L at constant density and find the point where 
they cross each other. Then, we use those points to con- 
struct the phase diagram. 



III. RESULTS 

A. Collective motion 

We will illustrate the collective dynamics in our model 
by sequentially changing one of the two main parameters: 
the density p and the parameter controlling energy influx 
rate, q, which therefore determines the average propul- 
sion speed of the particles. In Fig. [5] we display simula- 
tion snapshots obtained at a fixed input power q ~ 0.3 
and different particle number densities. All of the snap- 
shots feature significant density fluctuations. There is 
no obvious global ordering in the system but we can 
detect an increase of the cluster size. At fixed power 
q, an increase of the particle number density p leads to 
stronger density fluctuations and the velocity alignment. 
In Fig. [llb),(c), we can notice a formation of dense parti- 
cle groups at high concentrations, which move roughly in 
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the same direction. At low density (Fig. [IJa)), however, 
particles cannot form large groups and their velocities 
distributed randomly. 

Figure [3] illustrates the variation of the statistical char- 
acteristics of the swarm upon a change of the particle con- 
centration. The velocity autocorrelation function, C{t), 
in Fig. [3lja) shows an exponential decay at low density, 
p < 0.2, with the decay time increasing with the concen- 
tration. At the higher densities, the decay changes dra- 
matically, so the particles' direction of motion is getting 
much more stable in time. The spatial velocity correla- 
tion function, C\\ (r), in Fig. ^h) shows two distinct types 
of behavior: the decay is exponential at the two lowest 
concentrations, p = 0.1 and 0.2, while it becomes alge- 
braic at p > 0.2. We previously observed the transition 
to the power law form for two-point velocity correlations 
for the Vicsek-type model f39| . 

The cluster statistics for the density series is shown in 
Fig. m The plot in the inset confirms our observation 
that the cluster size is growing fast with the concentra- 
tion. At p = 0.05 and 0.1 the cluster size distribution, 
as shown in the main plot, decays exponentially while at 
larger density p it changes into a power law, which has 
been observed previously and is characteristic for the or- 
dered phase di-Iill. 

Now, we will look at the behavior of the system at 
constant density p = 0.2 while varying the energy in- 
flux rate q. As the ABP system changes the behavior 
from dissipative to driven upon increase of the energy 
influx rate, we expect the disordered motion at low q 
and onset of ordered behavior at high q levels. Fig. [5] 
shows the alignment of particles velocities at different in- 
put powers of the motor. At low q {q = 0.2, Fig.[5][a)) we 
see a practically homogeneous disordered system. Then, 
at g = 0.5 (Fig. mb)) well-shaped clusters are formed, 
within which the particles move in nearly the same di- 
rection. At high q, Fig. [51Jc), the velocities are high, the 
clusters are compact, and we observe a significant degree 
of alignment. Note also that in the system with g = 10 
there are very few single particles and most particles be- 
long to the swarm. Note also the shape of the swarm: the 
ABP with dissipative interactions tend to form bands, 
which are perpendicular to their velocity. 

Analysis of the velocity autocorrelation function in our 
model at fixed concentration p = 0.2 (Fig. Ula)) shows 
that at g > 0.3 the motion of particles is very persistent. 
The change of the direction of motion is realised only 
through collisions between different clusters. However, 
we see a fast decorrelation of the velocity at low q, q — 0.1 
and 0.3, due to the thermal noise. In the spatial velocity 
correlations (Fig. IH^b)) we see a sharp transition from 
the exponential decay at g = 0.1 and 0.3 to a power law 
decay for larger q. All the curves showing the power law 
decay are practically identical. 

Figure [7] presents the cluster statistics for the same 
density p = 0.2. The main plot shows the cluster size 
distribution on a log-log scale. We observe two qualita- 
tively different distributions: at low energy influx rates. 



FIG. 2: Typical distribution of particles inside a simulation 
box at constant propulsion power q = 0.3 and different parti- 
cle number densities: (a) p — 0.1, (b) p = 0.3, (c) p = 0.5. 
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FIG. 3: Statistical properties of the ABP-DPD model at a 
constant energy influx rate q = 0.3. (a) Semi-log plot of 
velocity autocorrelation function C{t) over time t. (b) Spatial 
velocity correlation function C||(r). 



q = 0.1 to 0.3; the curves show an exponential decay. At 
the higher q all of them are practically identical and have 
a straight segment at large numbers, which indicates the 
power law decay of the distribution. The transition can 
be located on the inset, where the evolution of the mean 
cluster size is shown. We see a kink on the curve at 
0.4. 

We measured the decay exponent for the spatial ve- 
locity correlation function, C||(r), in the whole range of 
studied parameters p and q and found that the exponent 
assumes universal values that depend only on the density 
but not on q, friction coefBcients or temperatures. The 
values of the exponent are plotted in Fig. |8l 



B. Orientational ordering 

Behaviour of the order parameter at various densities 
is plotted on Fig. At low densities the order param- 
eter values are close to zero during the whole simulation 
time which means that particle velocities are globally dis- 
aligned. At small propulsive power, q = 0.3, the ordering 
sets in slowly and reaches high values of about (/) « 0.8 
only at overlap densities of p w 1. At very high driving 
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FIG. 4: Cluster statistics for the ABP-DPD model at a con- 
stant energy influx rate q = 0.3. Inset: The average cluster 
size. The exponent pirn) oc m"'' for the straight segment: 
p = 0.3 C ~ 0.94, p = 0.4 C ~1, /5 = 0.5 C ~ IT, p = 0.6 
C^l.3. 



power, q = 100, the order parameter reaches unity at 
densities of about p = 0.04 that corresponds to the mean 
distance between the particles is r = p^^/^ = 5, which 
is greater than the radius of interaction ~ 3. Obvi- 
ously, the cohesive effect of the collisions keeps particles 
together, as can be seen already from the snapshots in 
Fig. El 

The phase diagram for our system is shown on Fig. 
[TOl The location of the transition points for each set of 
parameters was determined using the standard Binder 
cumulant analysis from the intersection of the cumulant 
curves Gl calculated for three different system sizes. The 
behavior of the cumulant indicates the continuous char- 
acter of the transition. It is clearly seen that the ordered 
behavior, at fixed environmental noise, is possible at cer- 
tain minimum energy influx rate g, which, in its turn, 
determines the average propulsion speed. The critical 
energy influx rate changes with concentration according 
to the power law Qc oc p~'^ , where k is 0.46 ± 0.02. We 
show the transition lines for two ambient temperatures, 
= 0.3 and 0.6. The twice as higher temperature of 
the environment at fixed friction 7^ means that the pas- 
sive fluctuations {D^) are twice as more intense and a 
higher energy influx is required for the ABP to be able 
to align. The qc values required for the transition at 
~ 0.6 and = are roughly 1.4 times higher than 
those found at = 0.3 and T'^ = 0. The q^ values 
observed =T^ = 0.3 are very close to those obtained 
at = 0.6 and = 0. 

Here, we would also like to demonstrate how the par- 
ticle individual and collective dynamics depends on the 
key parameters of the interaction. As mentioned above, 
the velocity correlations in our system decay exponen- 
tially in time in both phases according to C{t) oc e""*/"^ 
(see Fig. [TT|) . We have measured the correlation time 
to demonstrate the role of the intraswarm dissipation, 
which is controlled by 7'^. In a system without inter- 
actions, the relaxation time would be completely deter- 
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FIG. 5: Typical distribution of particles inside simulation box 
at constant density p — 0.2. (a) q = 0.2, (b) g — 0.5, (c) 
q = 10. 
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FIG. 6: Statistical properties of the ABP-DPD model at con- 
stant density p = 0.2. (a) Semi-log plot of velocity auto- 
correlation function C(t) over time t. (b) Spatial velocity 
correlation function C\\(t\ 




FIG. 7: Cluster statistics for ABP-DPD model at constant 
density (p = 0.2). Inset: Average cluster size. The exponent 
p{m) cx m""- for the straight segment: q = 0.7 C, ~ 1.16, 
g = 1.0 C 1-1, g = 5.0 C « 1-1, q = 10.0 C ~ l-ll- 



mined by the dissipative and driving mechanisms of the 
Langevin equation and would normally decrease with in- 
creasing the friction, = 771/7^. In contrast, as can be 
seen in the plot, the correlation time in the swarm, r'^, is 
growing proportionally to 7"^. The friction coefficient 7'^ 
scales the dissipative power of the pairwisc collisions and 
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FIG. 8: Behavior of the exponent for the velocity correlation 
function, C||(r) oc r-'*+=-'' 



FIG. 10: Phase diagram for the ABP-DPD model at 7^ = 0.3, 
7^ = 1, = 0. 




FIG. 9: Orientational order parameter for our model {ip = 1 
corresponds to completely ordered system, </p = - to com- 
pletely disordered system). 



therefore is the key parameter controlhng the alignment. 
Moreover, in the main plot wc sec that the mean order 
parameter is also growing larger with 7"^. 

Finally, Fig. [12] illustrates the role of the swarm tem- 
perature . The swarm temperature in our model can 
be defined via a fluctuation-dissipation relation for the 
parameters of the pairwise interaction, noise and the fric- 




FIG. 11: Orientational order parameter in the ABP-DPD 
model as a function of intraswarm friction 7"^. Inset: Velocity 
correlation time vs friction 7'^ {T^ — 0, p — 0.2, q — 0.3). 



tion coefficient, as given by the Eq. ([TT|) . In terms of the 
temperature, the transition looks completely analogous 
to what is seen in the magnetic systems. At zero tem- 
perature, the ordering is maximal, while it is suppressed 
by the fluctuations and vanishes at certain maximal tem- 
perature , which is increasing with the input power 
q. The order parameter approaches zero according to a 
power law oc \T^ - T^f with /^=0.52 for q = 0.5, 
/3 = 0.41 for q ~ 0.7, and (3 ~ 0.37 for q ~ I, which is 
in agreement with the critical exponent /S reported ear- 
lier for the Vicsek model and other models with aligning 
interactions Ell HI. 



IV. DISCUSSION 

As we can see from the numerical data, a system of 
ABP with dissipative interactions indeed demonstrates 
the same qualitative properties as the well studied Vic- 
sek model. In contrast to most previous approaches, 
here we have separated two influences on the particle 
motion: the effect of the environment, which is intro- 
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FIG. 12: Behaviour of the order parameter as a function of 
the swarm temperature (7^ = 1, P = 0.2, = 0.3). The 
dashed lines show the power law fit to the points left of the 
transition temperatures. 



duced through the Langcvin equation ([T]), and the effect 
of the other active agents, where the interactions are set 
by its own pairwise dissipative parameter. Both effects 
can be associated with a temperature, a corresponding 
noise and a friction parameter that controls the rate of 
dissipation. Note that these two types of noise and dissi- 
pation have different influence on the system. The former 
one is acting even on single particles, while the latter ap- 
plies only to the groups and automatically vanishes for 
single agents. At zero temperature of the environment, 
= 0, the model reduces to the motion with a constant 
speed, usually referred to as self-propelled particles, as 
for instance in the Vicsek model. At = 0, we have a 
system with aligning inelastic collisions but without the 
corresponding active noise. As we see from Fig. [Til the 
pairwise friction that scales the dissipation power in the 
collisions can be used to regulate the amount of alignment 
in the system, and, hence, the mean order parameter. 
Thus, our model allows one to mix these contributions in 
different proportions and model different swarming sce- 
narios. 

Now, we would like to discuss the extent of the dif- 
ferences and similarities of the swarming behavior in our 
model to observations from the Vicsek model in more de- 
tail. The main difference of the present analysis from the 
previous studies is that in most simulations we assumed 
a constant noise, as it is associated with the action of the 
environment, and varied the propulsive power of the par- 
ticles. This path, however, can be conveniently mapped 
onto a situation with a fixed particle speed and a variable 
noise. In static isolated systems, the ratio of the char- 
acteristic interaction energy to the thermal energy com- 
pletely determines the equilibrium state. In the swarm 
of active particles, the crucial number is the ratio of the 
stationary velocity due to propulsion to the characteris- 
tic velocity due to thermal fluctuations. This ratio can 



1/2 



(20) 



where we used the relation between the friction and the 
relaxation time in the Langevin equation, = 
and the equipartition relation, m{V^q) = ksT^ , Thus, 
the ratio in question is equivalent to the incoming en- 
ergy within the characteristic relaxation time, qr^ ^ to 
the thermal energy. In case the noise level is fixed by 
T^, is it the stationary particle speed that matters. The 
mean speed in the ABP model at large q is given by 
~ 9/7^- At higher temperatures of the environment, 
one needs to pump in more energy to produce the same 
ratio of the characteristic speeds. This point is confirmed 
by the data presented in Fig. (TU] In case the station- 
ary speed is fixed, one needs to reduce the temperature, 
which is, in the Langevin or DPD models, proportional 
to the fluctuation amplitude D. Therefore, the phase di- 
agram in terms of noise amplitude D vs density p or 
vs p will be inverse of our diagram shown in Fig. 1101 It 
is interesting to note that the sum of the critical swarm 
temperature, as shown in Fig. I12[ and the ambient tem- 
perature is roughly proportional to the energy influx 
rate. Here, we have ^T^ w 0.2+0.3 = 0.5 for q = 0.5, 
T^ + T^ « 0.7 for q ^ 0.7, and + w 0.8 + 0.3 = 1.1 
for q — 1. So, the fluctuations of different nature sim- 
ply add up to increase the effective swarm's temperature, 
which can be defined as T = + . This idea is sup- 
ported also by the phase diagrams shown in Fig. [TUJ 
where two systems with equal values of + = 0.6 
demonstrate a transition at nearly the same q and p. 

In simulations, we observe aggregation and orienta- 
tional ordering of ABP at sufficiently high densities in 
presence of sufficiently high propulsive power. At the 
fixed level of noise and propulsive power, the cluster size 
grows with the particle concentration in the same way as 
we observed previously for the Vicsek model [s^. Thus, 
the dissipative interactions as well lead to cohesion of ac- 
tive particles. Secondly, they lead to particle alignment 
as can be seen from the growth of the order parameter 
with concentration, again, similar to the dependence seen 
in the Vicsek model. The transition into the orientation- 
ally ordered phase happens along the line qc oc 
which is an inverse of the transition line for the Vicsek 
model, where it happens at ^ oc where ^ is the noise 

amplitude [H, . This power law behavior seems to be 
not unique to the Vicsek model, but a universal property 
of systems with global alignment and has been reported 
also for systems with pairwise alig ning interactions (with 
an exponent k = 0.46 ± 0.04) [26f. 

In what regards other properties, we should mention 
the behavior of the correlations functions C {t) and C|| (r) 
(Figs, inland [3]), which demonstrate the same qualitative 
features as the Vicsek model we studied previously [39| . 
The two-point velocity correlation function changes the 
shape from exponential to a power law at the critical 
point and inside the whole region of the ordered behavior. 
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The exponent rj, which describes the decay of C||(r) oc 
-j-jjg ordered phase, takes the same values from 
0.5 to 0.97 on increasing density and shows the same 
density dependence as we previously saw in the Vicsek- 
type model. It seems to be insensitive to other details of 
the system and reflects just the nature of the ordering. At 
the transition point, the exponent is expected to satisfy 
the Fisher's scaling law: 7/1^ = 2 — rj [44|, where 7 and 
v are the critical exponents for isothermal succeptibility 
and the fluctuation correlation radius. In the limit of low 
concentrations, where the repulsions are not important, 
we have 77 = 0.5, thus 2 — r/ = 1.5, which is in agreement 
with the result for 7/1/ = 1.47 obtained previously for 
the standard 2D Vicsek model [42j. Moreover, the shape 
of the cluster size distributions as shown in Figs. [7] and 
13] in our model is also identical to that for the Vicsek 
model ranging from ( ~ 0.5 to C = 1-5 depending on the 
level of noise and the density [39l - l41| . Although the type 
of active particle, the interactions and the type of noise 
differ from those in the Vicsek model, the identical values 
of the exponents suggest that our model belongs to the 
same universality class [4l| . 



particle interactions, which are introduced via inelastic 
collisions. We find that the ABP-DPD model exhibits 
an orientational order-disorder transition on increasing 
energy influx rate or particle number density, which is 
completely analogous to that in the Vicsek model. We 
have shown that the parameter space of such an active 
system can be reformulated in terms of effective temper- 
atures of the environment and of the swarm and the ratio 
of the characteristic thermal energy to the energy influx 
per particle. 
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V. CONCLUSIONS 



We have studied dynamic self-organization in a model 
combining the active Brownian particles with dissipative 
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